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Executive  Summary 

The  overall  goal  of  this  research  is  to  contribute  to  national  defense  and  security  by  advancing 
the  state-of-the-art  in  explosive  dispersal  of  particles/droplets  and  replacing  current  empiricism 
with  a  physics-based  fundamental  approach.  This  will  enable  game-changing  advancements  in  (i) 
lethality  technologies  for  the  development  of  next  generation  tunable-output  explosives,  (ii) 
force-protection  technologies  for  defeating  improvised  explosive  devices  and  mines,  and  (iii) 
counter-terrorism  technologies  for  neutralizing  chemical  and  biological  weapons  of  mass 
destruction. 

The  core  intellectual  objective  is  to  radically  advance  the  field  of  Compressible  Multiphase 
Turbulence  (CMT)  by  systematically  answering  key  outstanding  questions.  We  will  advance  a 
rigorous  physics-based  fonnulation  of  CMT  and  thereby  enable  accurate  high-fidelity 
simulations  with  true  predictive  capability.  The  longer-term  technological  objective  is  to 
facilitate  the  development  of  a  Virtual  Explosive  Dispersal  Facility  (VEDF),  which  can  be  used 
to  perform  physics-based  simulations  of  explosive  dispersal  of  particles/droplets  in  applications 
of  interest  to  DoD. 

We  have  focused  on  the  following  two  key  questions:  (1)  Can  the  current  modeling  and 
simulation  approaches  accurately  predict  explosive  dispersal  particles  observed  in  the  above 
examples  (and  other  related  experiments)?  (2)  If  not,  what  additional  fundamental  understanding 
is  needed  that  can  be  translated  to  improved  modeling  and  simulation  capabilities? 

Existing  simulation  tools  have  become  quite  advanced  over  the  past  decade  and  they  incorporate 
a  wide  variety  of  capabilities.  Despite  these  advancements,  their  predictive  capability,  especially 
for  the  class  of  problems  of  present  interest,  has  been  quite  poor.  The  fundamental  difficulty  is  in 
lack  of  complete  understanding  of  the  complex  physics  of  multiphase  flow  under  conditions  of 
extreme  compressible  turbulence.  In  order  to  accomplish  these  intellectual  and  technological 
objectives,  we  have  studied  in  detail  the  following  micro  and  mesoscale  problems: 

•  Investigate  at  the  microscale  the  complex  thermo-mechanical  interaction  between  an  isolated 
particle  and  the  surrounding  rapidly  expanding  compressible  flow.  Develop  rigorous  physics- 
based  closures  for  interphase  coupling  between  the  particle  and  the  gas. 

•  Extend  the  interphase  coupling  model  to  finite  Reynolds  and  Mach  numbers,  and  to  account 
for  the  effects  of  deformation  and  compressive  heating  in  case  of  strong-shock-particle 
interaction.  Our  initial  effort  was  to  develop  and  test  a  point-particle  model  to  accurately 
predict  the  mass-average  velocity  and  temperature  of  the  particle.  However,  from  the  direct 
numerical  simulations  it  was  noted  that  the  peak  values  of  stress  and  temperature  are  far 
higher  than  the  average  values.  Efforts  are  underway  to  extend  the  point-particle  models  to 
predict  these  peak  values  as  well. 

•  Extend  the  interphase  coupling  models  from  an  isolated  particle  to  a  cluster  of  particles  by 
incorporating  the  effect  of  finite  volume  fraction.  Here  both  the  limits  of  dilute  and  dense 
volume  fractions  must  be  considered.  In  the  former  only  the  indirect  particle-particle 
interaction  as  mediated  by  the  surrounding  gas  will  be  important,  while  in  the  later  limit 
additional  effects  of  direct  particle-particle  collisional  interaction  must  be  modeled.  This  work 
is  in  collaboration  with  Drs.  J.L.  Wagner,  S.J.  Beresh,  and  S.P.  Kearney  at  the  Sandia 
National  Laboratories. 

•  Existing  compressible  multiphase  flow  formulations  are  not  universal.  There  are  two  classes 
of  formulations  in  use.  The  first  class  of  fonnulations  apply  to  situations  of  dense 
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compressible  multiphase  flows,  where  compression  of  the  dispersed  phase  plays  an  important 
role.  Multiphase  formulations  of  this  type  follow  along  the  lines  of  the  original  work  by  Baer 
&  Nunziato  [1].  The  later  class  of  formulations  are  appropriate  for  less  dense  and  less  intense 
scenarios,  where  the  compressional  effect  on  the  particles  can  be  ignored.  Here  we  seek  to 
develop  a  unified  multiphase  formulation  that  seamlessly  handless  both  these  limits 
accurately.  Such  a  formulation  is  of  critical  importance,  since  explosive  dispersal  problems  of 
present  interest  start  from  a  dense  flow  limit  to  quickly  evolve  into  a  dilute  particulate  flow. 
This  work  has  been  initiated  in  collaboration  with  Dr.  T.P.  McGrath  at  NSWC-Indian  Head, 
MD. 

•  Rigorously  apply  the  interphase  momentum  and  energy  coupling  models  to  problems  of 
shock-particle  and  detonation-particle  interaction  and  validate  the  coupling  model  against 
available  high-quality  experimental  results.  In  particular,  validation  of  thermal  coupling 
models  is  of  great  interest,  as  there  are  not  adequate  validation-quality  experimental 
measurements  of  particle  and  gas  temperature  immediately  following  a  strong-shock  or 
detonation  particle  interaction. 

•  Investigate  the  process  of  explosive  dispersal  of  particles,  starting  from  the  initial  instability 
of  a  rapidly  expanding  layer  to  its  subsequent  dispersal  process.  Here  the  initial  inviscid 
instability  analysis  was  first  extended  to  viscous  instability  analysis,  where  the  effect  of 
viscosity  is  to  identify  a  most  amplified  spherical  mode.  We  are  in  the  process  of  extending 
the  gas-gas  contact  interface  instability  analysis  to  multiphase  particulate  fronts.  An  important 
advancement  in  the  stability  analysis  is  the  realization  of  the  key  role  of  interface  thickness  in 
dictating  the  wave  number  of  the  most  amplified  mode.  The  linear  stability  analysis  has  also 
been  extended  to  nonlinear  stages  of  bubble  and  spike  growth  regime  and  investigated  through 
direct  numerical  simulations. 

•  Rigorously  validate  the  micro  and  macroscale  formulations  and  results  against  high-quality 
experimental  data  available  at  the  Air  Force  Research  Laboratory  Munitions  Directorate  and 
at  the  DOE  NNSA  Laboratories. 

•  Train  students  and  postdocs  in  problem  of  defense  and  security  science  and  engineering. 
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2.  Introduction 

Lethality  technologies  to  develop  next  generation  tunable-output  explosives,  force-protection 
technologies  to  defeat  improvised  explosive  devices  (IEDs)Anines,  and  counter-terrorism 
technologies  to  neutralize  chemical  and  biological  weapons  of  mass  destruction  (WMD)  are  of 
great  importance  to  our  Nation's  security  and  defense.  Fundamental  to  these  technologies  are  the 
ability  to  control  explosive  dispersal  of  particulates/droplets,  to  direct  the  location,  time  and  rate 
of  energy  release,  and  to  mitigate  the  effects  of  blast,  fireball  and  impacting  projectiles. 

It  is  now  possible,  and  will  likely  become  increasingly  so  in  the  future,  to  produce 
multifunctional  explosives  whose  microstructural  details  of  particle  size  distribution  and 
chemical  composition  can  be  controlled  to  great  precision.  Thus,  there  is  clear  opportunity  for 
engineering  next  generation  tunable  weapons  systems  that  can  deliver  a  controlled  amount  of 
energy,  to  an  exact  location,  at  a  precise  moment,  and  at  a  desired  rate.  Tunable  explosives, 
however,  have  not  yet  been  fully  realized  in  practical  systems.  The  key  missing  link  is  to  connect 
the  microstructural  details  of  the  explosive  to  the  resulting  spatio-temporal  details  of  macroscale 
dispersion,  turbulent  mixing  and  volumetric  energy  release.  If  this  challenge  can  be  addressed 
game-changing  capabilities,  such  as  quantum  increase  in  energy  release  or  localized  release  for 
low  collateral  damage,  can  be  realized. 

Compressible  multiphase  turbulence  is  at  the  heart  of  this  problem.  Advances  in  our  fundamental 
understanding  of  compressible  turbulent  multiphase  flows  will  positively  impact  many  other 
problems  of  interest  to  DoD.  For  example,  the  devastating  effects  of  IEDs  arise  from  the 
fragments,  debris  and  incendiaries  that  they  explosively  disperse.  To  effectively  protect  against 
these  projectiles  it  is  important  to  accurately  know  the  momentum  and  energy  imparted  to  them 
in  an  explosion.  Or  in  the  case  of  blast  mitigation  with  water  mist,  it  is  the  complex  thermo¬ 
mechanical  interaction  of  the  water  droplets  with  the  compressible  turbulent  flow  that  mitigates 
the  blast. 

Compressible  multiphase  turbulence  is  in  its  infancy  and  poses  a  grand  challenge  as  it  combines 
three  complex  physics:  turbulence,  compressibility  and  multiphase  flow.  The  urgent  need  for 
better  fundamental  understanding  is  clear  from  the  national  defense  and  security  applications. 
Although  substantial  progress  has  occurred  over  the  past  several  decades  in  the  individual  areas 
of  turbulence,  compressible  and  multiphase  flows,  our  understanding  when  all  three  mechanisms 
are  simultaneously  active  is  rudimentary. 

The  primary  objective  of  the  ongoing  work  is  to  radically  advance  the  field  of  compressible 
multiphase  turbulence  by  systematically  answering  key  outstanding  questions.  We  will  advance  a 
rigorous  physics-based  formulation  and  thereby  enable  accurate  high-fidelity  simulations  with 
true  predictive  capability.  High  quality  experiments  that  have  been  performed  (or  being 
performed)  at  DoD  and  DOE  will  be  used  to  validate  the  models  and  the  simulation  framework. 

In  particular,  we  will  focus  on  the  class  of  problems  involving  explosive  dispersal  of  particles 
and  droplets.  These  dispersed  particles  can  be  inert  as  in  the  case  of  casing  fragments  or  liner 
elements  or  can  be  reactive  as  in  the  case  of  multiphase  explosives.  Using  fundamental  advances 
in  modeling  and  simulation  of  compressible  multiphase  turbulence  we  perform  predictive 
physics-based  simulations  of  selected  explosive  dispersal  problems.  By  comparing  results  with 
experimental  observations  and  measurements  the  predictive  capability  of  the  new  formulation  is 
evaluated. 
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3.  Scope  and  Challenge 

Consider  the  following  few  examples  of  explosive  dispersal.  In  Figure  1  high  speed  images  from 
a  thermobaric  test  is  shown.  In  this  case  there  are  no  particles  being  dispersed  and  therefore  can 
be  considered  as  single-phase  baseline.  The  instability  and  turbulent  nature  of  the  expanding  gas 
front  can  be  clearly  observed.  These  instabilities  are  thought  to  be  associated  with  Rayleigh- 
Taylor  (RT)  and  Richtmeyer-Meshov  (RM)  instabilities.  But  the  true  nature  of  the  instability  and 
transition  to  turbulence  at  such  rapidly  expanding  material  fronts  is  poorly  understood. 

Figure  2  shows  images  of  an  expanding  cloud  of  aluminum  particles  (Particle  image  velocimetry 
of  this  field  were  also  taken).  In  these  experiments,  performed  at  ARFL  Munitions  Directorate  at 
Eglin  Air  Force  Base,  four  stacked  A- 5  pellets  (see  frame  e)  were  surrounded  by  a  paper  jacket 
filled  with  aluminum  particles  and  the  resulting  pattern  of  explosive  dispersal  of  particles  was 
observed  and  measured.  These  experiments  were  repeated  for  aluminum  and  tungsten  particles 
of  varying  particle  size  (see  [2]  for  additional  details). 


Fig  1:  High  speed  image  of  a  thermobaric  test  of  3kg  TBX.  The  highly  turbulent  nature  of  the  expanding 
gas  contact  and  the  associated  fireball  can  be  clearly  seen  (from  Dr.  M.R.  Baer). 


Fig  2:  Early  time  expansion  (a)  An  H-10  powder  annulus  charge  at  50ps  after  initiation  from  the  right 
side,  (one  inch  per  block)  (b).  95  powder  annulus  charge  at  5 Ops  after  initiation  from  the  right  (c) 
tungsten  particle  charge  at  50ps  into  expansion,  no  combustion  due  to  the  inert  metal  powder  (d)  HSFC 
speed  images  of  TSS  case  breakup  with  fragment  images  at  50  ps  post  detonation  (from  [2]). 

Figure  3  shows  an  explosive  spherical  dispersal  of  glass  particles.  Prof.  Frost  and  co-workers 
repeated  these  experiments  under  a  range  of  conditions  involving  different  kinds  of  particles  and 
geometric  configurations.  In  particular,  we  draw  attention  to  the  fact  that  even  the  qualitative 
behavior  of  instability  of  the  material  front  and  the  pattern  of  particle  dispersion  in  these  three 
examples  are  quite  different.  Furthermore,  a  detailed  look  at  these  experimental  results  shows 
substantial  quantitative  variation,  which  cannot  be  explained  with  current  understanding. 
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Fig  3:  Spherical  dispersal  of  glass  particles  with  C4  explosive  using  120  mm  glass  particles  and  28g  of 
C4.  Note  the  needle  shaped  instability  of  the  expanding  particulate  front,  (from  [3]). 


We  now  ask  the  following  key  questions:  (1)  Can  the  current  modeling  and  simulation 
approaches  accurately  predict  explosive  dispersal  particles  observed  in  the  above  examples 
(and  other  related  experiments)?  (2)  If  not,  what  additional  fundamental  understanding  is 
needed  that  can  be  translated  to  improved  modeling  and  simulation  capabilities? 


Existing  simulation  tools,  both  commercial  as  well  as  those  developed  and  used  by  individual 
research  groups,  have  become  quite  advanced  over  the  past  decade  and  they  typically  incorporate 
a  wide  variety  of  capabilities  to  handle  complex  geometries,  particulates  and  droplets,  ambient 
flow  turbulence  (through  Reynolds-averaged  Navier-Stokes  and  large-eddy  simulation 
techniques).  Despite  these  capabilities,  their  predictive  quality,  especially  for  the  class  of 
problems  of  present  interest,  has  been  quite  poor.  The  fundamental  difficulty  is  in  lack  of 
complete  understanding  of  the  complex  physics  of  multiphase  flow  under  conditions  of  extreme 
compressible  turbulence.  The  coupling  models  that  are  in  current  use  for  mass,  momentum  and 
energy  transfer  between  the  particles/droplets  and  the  surrounding  gas  are  inadequate,  as  they  are 
generally  based  on  our  understanding  and  modeling  for  simple  steady,  low-speed  conditions. 

The  problem  of  explosive  dispersal  of  particles  requires  us  to 

•  Develop  a  universal  set  of  governing  equations  for  compressible  turbulent  multiphase  flows 
that  is  applicable  over  the  entire  range  of  particulate  concentration  from  the  dilute  limit  to 
very  dense  suspension  and  from  low-speed  incompressible  flows  to  very  high-speed  high 
density  flows. 

•  Carefully  address  interphase  coupling  at  conditions  that  are  of  relevance  to  DoD  applications 
and  to  develop  physics-based  improved  models  of  momentum  and  energy  coupling.  These 
coupling  models  will  appear  as  closure  relations  in  the  universal  set  of  governing  equations 
addressed  above.  The  interphase  interaction  is  significantly  complicated  due  to  (a)  the  highly 
unsteady  nature  of  the  problem,  (b)  interaction  of  compressible  flow  features  such  as  shocks 
and  contacts  with  the  particles,  (c)  high  Reynolds  number,  (d)  compressibility  effects  of  high 
Mach  number,  and  (e)  particle-particle  interaction  at  large  volume  fractions. 

•  Address  additional  complications  that  arise  from  random  particle  size,  shape  and  distribution. 

•  Investigate  the  instability  of  rapidly  expanding  gas  and  particulate  fronts  and  study  the  nature 
of  flow  transition  and  development  of  multiphase  turbulence. 
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4.  Results  and  Discussion 


The  results  of  the  present  work  will  be  presented  under  five  different  topics.  Our  efforts  in  all 
these  five  topics  have  produced  results  that  have  already  appeared  as  journal  articles  (or  in 
review  process  or  being  prepared).  These  archival  journal  publications  are  listed  as  items  J1  to 
J16  in  section  6.  Therefore  the  results  and  discussions  to  be  presented  here  will  be  brief,  with 
emphasis  on  only  key  sample  results.  The  details  can  be  obtained  from  the  journal  articles.  In 
section  5  we  will  discuss  ongoing  and  future  research  efforts  that  closely  follow  those  listed  here. 

4.1  Theoretical  Foundation  -  Exact  Solutions  For  an  Isolated  Particle 

Consider  the  scenarios  shown  in  Figures  2  and  3.  In  each  case,  the  dispersive  pattern  of  millions 
of  particles  is  desired.  At  this  macroscale  it  will  be  impossible  to  resolve  the  flow  around  each 
and  every  particle.  One  needs  to  resort  to  a  point-particle  approach,  where  each  particle  is  treated 
as  point  and  its  motion  is  computed  based  on  aerodynamic  forces  that  act  on  it.  In  the  context  of 
extreme  compressible  flows  under  consideration,  the  force  on  the  particle  will  be  a  function  of 
(at  the  least)  (1)  particle  velocity,  (2)  particle  acceleration,  (3)  local  fluid  velocity,  (4)  local  fluid 
acceleration,  (5)  local  particle  volume  fraction,  and  (5)  other  local  fluid  properties  such  as 
density,  pressure  gradient  and  speed  of  sound. 

While  at  low-speeds  (in  the  incompressible  regime)  one  can  correlate  drag  coefficient  as  a 
function  of  only  Reynolds  number  based  on  relative  velocity,  this  simple  approach  is  grossly 
inadequate  at  the  extreme  conditions  of  interest  here.  All  the  other  parameters,  especially  effects 
of  unsteadiness,  compressibility  and  volume  fraction  will  play  an  important  role.  However,  a 
naive  attempt  at  correlating  the  force  on  the  particle  in  terms  of  the  long  list  of  parameters  will  be 
futile. 

A  physics-based  approach  to  modeling  this  interphase  momentum  and  energy  coupling  is 
essential  for  success.  As  a  theoretical  foundation,  we  take  the  limit  of  zero  particle  Reynolds  and 
Mach  numbers.  In  this  limit  the  linearized  compressible  Navier-Stokes  equations  of  flow  around 
a  particle  can  be  rigorously  solved.  The  following  steps  have  allowed  us  to  establish  foundational 
results  which  will  then  be  exploited  in  our  other  research: 

•  The  problem  of  an  isolated  sphere  undergoing  arbitrary  motion  in  a  quiescent  compressible 
fluid  is  first  solved  to  obtain  an  exact  expression  for  force.  Individual  mechanisms  in  the  form 
of  quasi-steady,  inviscid-unsteady,  and  viscous-unsteady  contributions  to  overall  force  are 
identified.  Numerical  results  are  in  perfect  agreement  with  the  theory  (see  Figure  4).  The 
resulting  interphase  coupling  fonnulation  offers  an  explicit  expression  for  the  unsteady  force 
in  the  time  domain  and  can  be  considered  as  a  generalization  of  the  classic  Basset- 
Boussinesq-Oseen  (BBO)  equation  to  compressible  flow  (See  journal  article  J1  for  details). 

•  The  above  theoretical  analysis  was  extended  to  the  more  general  case  of  an  isolated  sphere 
undergoing  arbitrary  motion  in  a  spatially-varying  time-dependent  compressible  ambient 
flow.  The  problem  was  formulated  in  a  reference  frame  attached  to  the  particle  and  force 
contribution  from  the  undisturbed  ambient  flow  and  the  perturbation  flow  were  separated. 
Using  density  weighted  velocity  transformation  and  reciprocal  relation,  the  total  force  was 
first  obtained  in  the  Laplace  domain  and  then  transformed  to  the  time  domain.  Total  force  was 
again  separated  into  the  quasi-steady,  inviscid-unsteady,  and  viscous-unsteady  contributions. 
The  above  rigorously  derived  particle  equation  of  motion  can  be  considered  as  the 
compressible  extension  of  the  classic  Maxey-Riley-Gatignol  (MRG)  equation  of  motion.  It 


PI:  S.  Balachandar,  University  of  Florida 


7  |  P  a  g  e 


rigorously  incorporates  important  physics  that  arise  from  the  combined  effects  of 
inhomogeneity  and  compressibility  (See  journal  article  J5  for  details). 

•  Finally,  the  above  formulation  for  a  rigid  sphere  was  generalized  to  investigate  the 
corresponding  results  for  a  droplet.  Here  the  linearized  viscous  compressible  Navier-Stokes 
equations  are  solved  both  inside  and  outside  of  the  spherical  drop.  An  expression  of  the 
transient  force  is  first  obtained  in  the  Laplace  domain  and  then  transformed  to  the  time 
domain.  The  total  force  is  separated  into  the  quasi-steady,  the  inviscid-unsteady,  and  the 
viscous-unsteady  contributions.  The  new  force  expression  reduces  to  known  results  in  the 
limits  of  a  drop  in  an  incompressible  flow  or  a  rigid  particle  in  a  compressible  flow  (See 
journal  article  J7  for  details). 


Fig  4:  Time  evolution  of  the  normalized  unsteady  force.  Theoretical  predictions  are  plotted  as  solid  lines 
for  different  Mach  and  Reynolds  numbers.  Corresponding  simulation  results  for  the  different  cases  are 
shown  as  symbols.  Interestingly,  excellent  agreement  can  be  observed,  even  at  finite  values  of  Reynolds 
and  Mach  numbers.  Nonlinear  effect  can  be  shown  to  be  important  only  for  time  larger  than  (l/Re*M). 

•  It  is  well  known  that  the  computation  of  the  history  force  by  evaluating  the  convolution 
integral  is  very  demanding  in  terms  of  CPU  and  memory  requirements.  In  this  research  we 
seek  an  elegant  mathematical  approximation  that  greatly  simplifies  the  computation  of  the 
history  integrals  without  much  loss  of  accuracy.  This  approximation  will  greatly  benefit  the 
computational  accuracy  and  efficiency  of  Eulerian-Lagrangian  approach  to  multiphase  flows. 
Using  rational  approximation  for  the  history  kernel  in  the  frequency-space  we  have 
developed  an  approximate  formulation  of  the  viscous-unsteady  force  which  is 
computationally  far  more  efficient  than  the  traditional  convolution  integral.  For  the 
compressible  inviscid  unsteady  force,  the  kernel  in  the  frequency  domain  is  a  rational 
function.  As  a  result  the  convolution  integral  and  the  differential  formulation  of  the 
compressible  inviscid  unsteady  force  are  identical.  The  viscous-unsteady  kernel  in  the 
frequency-space  involves  fractional  powers  of  the  frequency  and  thus  only  an  approximate 
rational  form  can  be  obtained  using  N-tenn  Pade  approximations.  A  2-term  Pade 
approximation  has  been  shown  to  approximate  the  finite  Reynolds  number  viscous-unsteady 
kernel  to  within  1%  accuracy  for  all  non-dimensional  frequencies.  From  a  simple  analysis  of 
particle  motion  in  turbulent  flows,  we  argue  that  this  range  of  frequency  is  adequate  for  a 
wide  range  of  particle  properties  and  turbulence  Reynolds  numbers  (see  J1 1  for  further 
details). 
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4.2  Mesoscale  Mathematical  Formulation 

Recently,  in  collaboration  with  Dr.  T.P.  McGrath  at  the  Naval  Surface  Warfare  Center  -  Indian 
Head,  MD,  we  have  been  developing  a  rigorous  mathematical  formulation  for  the  compressible 
multiphase  flow  that  is  universal  in  its  applicability.  The  present  mathematical  formulations  of 
compressible  multiphase  flow  fall  within  one  of  two  categories.  The  first  kind  of  mathematical 
models  use  ideas  of  rational  mechanics  and  have  been  developed  for  dense  particulate  flows, 
where  the  dispersed  phase  undergoes  strong  compressional  effects.  But  the  compaction  equation 
for  particle  volume  fraction  employed  in  these  models  is  not  applicable  as  the  compression  effect 
becomes  negligible  as  the  volume  fraction  decreases.  The  second  kind  of  mathematical  models 
are  appropriate  for  dilute  particulate  flows,  but  cannot  account  for  compressional  effects  and  as  a 
result  fail  as  the  volume  fraction  increases  and  start  to  approach  packing  limit.  However,  for  the 
problem  of  present  interest  where  a  dense  packed  bed  of  particles  is  explosively  dispersed,  we 
require  a  universal  mathematical  formulation  that  spans  from  the  dense  regime  to  the  dilute 
regime. 

The  improved  universal  formulation  will  incorporate  the  following  important  features.  A  revised 
compaction  equation  as  shown  below  is  being  investigated 

^-  +  urVa  +  <l>V-ui  =  f(p-pd) 
at 

where  a  is  the  volume  fraction  of  the  particulate  phase,  u,  is  the  interface  velocity,  $  is  a 
weighting  factor,  whose  value  needs  to  be  chosen,  and  f(p~pd )  is  the  pressure  equilibration 
contribution  to  particle  compaction.  The  above  formulation  has  the  flexibility  to  reduce  to  the 
classic  expression  by  Baer  &  Nunziato  [1]  in  the  limit  ^  >  0  .  But  also  has  the  ability  to 
accommodate  constant  density  particulate  phase  in  the  limit  cj)  — >  a  .  We  are  in  the  process  of 
examining  the  thermodynamical  consequences  of  the  above  enhanced  compaction  equation. 
Also,  the  mathematical  structure  of  the  resulting  governing  equations  changes  and  the 
hyperbolicity  of  the  system  along  with  the  appropriate  definition  of  the  speed  of  sound  have  been 
derived.  Once  the  structure  of  the  universal  system  has  been  fully  established  we  will  perform 
one-dimensional  tests  and  then  proceed  to  implement  in  the  three-dimensional  code  of  our  group 
and  also  in  the  Navy  hydrocode  DYSMAS. 

4.3  Rigorous  Application  to  Shock-Particle  Interaction 

Here  we  build  upon  the  solid  foundation  laid  by  the  rigorous  theoretical  formulation  outlined  in 
section  4.1.  We  apply  the  interphase  coupling  formulation  to  the  problem  of  shock-particle 
interaction.  The  interaction  can  be  accurately  resolved  by  direct  numerical  simulations.  However, 
as  the  length  scales  of  interest  are  much  larger  than  the  particle  size  in  many  applications,  fully 
resolving  the  flow  around  millions  of  particles  is  impractical.  Therefore,  a  point-particle 
approach  is  inevitable  and  in  this  context  a  rigorous  model  for  momentum  and  energy  exchange 
is  very  important. 

Shock-particle  interaction  is  strongly  time-dependent,  so  unsteady  mechanisms  are  important  in 
momentum  and  energy  transfer.  The  theoretical  foundation  outlined  in  section  4.1  provides  us  a 
rigorous  physics-based  function  form  for  the  different  contributions  to  the  overall  force  and  heat 
transfer.  Here  we  extend  this  model  to  account  for  the  effects  of  finite  Reynolds  and  Mach 
numbers.  The  model  is  used  to  investigate  particle  interactions  with  a  planar  shock  wave  and  a 
spherical  shock  wave. 
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The  key  aspect  of  the  present  model  is  that  it  accurately  accounts  for  the  unsteady  mechanisms, 
which  are  often  ignored  in  the  traditional  models.  During  the  shot  duration  of  shock-particle 
interaction  the  fluid  acceleration  seen  by  the  particle  is  very  large  and  contributes  to  the  unsteady 
mechanisms.  The  peak  values  and  the  net  effects  of  unsteady  contributions  to  force  and  heat 
transfer  are  used  to  measure  the  importance  of  the  unsteady  contributions.  The  results  in  Figure  5 
show  the  ratio  of  peak  values  of  unsteady  contributions  of  force  (frame  a)  and  heat  transfer 
(frame  b)  to  their  quasi-steady  contributions  for  a  wide  range  of  particle  parameters.  It  is  clear 
that  for  a  wide  range  of  shock  Mach  numbers  and  particle  Reynolds  numbers  the  forces  and  heat 
transfer  due  to  unsteady  mechanisms  are  more  than  an  order  of  magnitude  larger  and  thus 
ignoring  them  can  result  in  very  large  errors.  Further  details  can  be  found  in  the  journal  article 
J2. 


Fig  5:  Here  we  consider  the  problem  of  force  and  heat  transfer  on  a  stationary  spherical  particle 
subjected  to  a  shock.  Frame  (a)  shows  contours  of  the  ratio  of peak  inviscid  unsteady  (pressure-gradient) 
force  to  the  peak  quasi-steady  force  (standard  drag  laws  are  limited  to  this  quasi-steady  force).  Frame  (b) 
shows  the  ratio  of  peak  inviscid  unsteady  heating  to  the  peak  quasi-steady  heat  transfer.  The  ratios  are 
plotted  for  a  range  of  shock  Mach  number  (x-axis)  and  for  a  range  of particle  size,  parameterized  as 
particle  Reynolds  number  (y-axis).  It  can  be  seen  that  the  ratio  can  be  as  large  as  20  for  the  case  of  force 
and  300 for  the  case  of  heating,  thus  demonstrating  the  importance  of  unsteady  mechanisms  over  a  wide 

range  of  shock-particle  interaction  conditions. 


Fig  6:  Here  we  consider  the  "spherical  multiphase  shock  tube  problem  ",  where  a  spherical  container  of 
unit  radius  is  initially  filled  with  a  high  pressure  gas-particle  mixture.  At  t  =  0  the  container  is  removed 
instantaneously  and  the  mixture  is  exposed  to  the  surroundings,  leading  to  a  blast  wave,  gas  contact, 
outward  propagating  front  of particles  and  an  expansion  fan  propagating  inward.  In  the  spherical 
geometry  the  gas  contact  reaches  an  asymptotic  final  radius,  which  is  presented  in  frame  (a).  Also  plotted 
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in  frame  (a)  is  the  location  of  the  particle  front,  as  a  function  of particle  diameter  (normalized  by  initial 
charge  radius),  at  three  different  times.  It  can  be  observed  that  (i)  particles  in  the  present  example  move 
past  the  gas  contact  into  the  shock-heated  ambient,  (ii)  There  is  an  optimal  particle  size  that  reaches  the 
farthest  at  any  given  time,  and  (Hi)  the  optimal  particle  size  increases  with  increasing  time.  In  Frame  (b) 
we  plot  the  error  in  particle  front  location  by  not  including  the  unsteady  mechanisms  for  varying  initial 
particle-to-fluid  density  ratio  (x-axis)  and  particle  size  (y-axis).  It  can  be  seen  that  the  error  in  particle 
front  location  can  be  as  large  as  5  initial  charge  radii,  and  from  frame  (a)  it  can  be  seen  that  this  error  is 
a  large  fraction  of  the  actual  gas  front  and  particle  front  location. 

We  extended  the  above  investigation  to  consider  the  interesting  phenomenon  of  particle  dispersal 
by  blast  waves.  A  model  problem,  i.e.,  a  sudden  release  of  a  compressed  gas-particle  mixture 
contained  in  a  spherical  container,  is  employed  to  investigate  the  fundamental  physics  of  particle 
dispersal.  The  problem  is  simulated  with  the  interphase  coupling  models  outlined  in  J2.  At  early 
times,  when  particles  are  accelerated  in  the  expansion  fan,  unsteady  force  and  heating 
contributions  are  much  larger  than  the  corresponding  quasi-steady  contributions.  Consequently, 
neglecting  unsteady  contributions  leads  to  significant  errors  in  particle  front  location  (the 
boundary  of  the  particle  cloud  -  see  Figure  6).  The  complex  wave  interactions  in  the  flow  have 
strong  influence  on  the  particle  motion.  As  a  result,  the  particle  motion  is  a  non-monotonic 
function  of  particle  density  or  diameter  and  the  evolution  of  particle  concentration  is  non- 
uniform  and  unsteady.  Further  details  can  be  found  in  the  journal  article  J3. 

4.4  Scaling  Analysis 

Modeling  of  inter-phase  coupling  is  critical  to  accurate  simulation  of  dispersed  multiphase  flows. 
In  this  effort,  we  address  three  important  fundamental  questions  about  inter-phase  coupling:  (1) 
Is  it  possible  to  simplify  the  integral  representation  of  viscous-unsteady  force?  (2)  Under  what 
conditions  unsteady  forces  are  important  in  forward  momentum  coupling?  (3)  What  parameters 
accurately  evaluate  the  back  effect  of  both  the  quasi-steady  and  unsteady  forces  in  the  fluid 
momentum  equation?  The  key  findings  of  the  present  study  in  regard  to  these  three  questions  are 
summarized  separately  as  follows.  The  details  of  the  analysis  can  be  found  in  J13-J15. 

Conclusions  on  the  representation  of  the  viscous-unsteady  force:  In  general,  the  viscous-unsteady 
force  is  given  by  a  Basset-like  convolution  integral,  whosecomputation  is  often  costly.  However, 
when  both  the  particle  and  the  surrounding  fluid  acceleration  are  sufficiently  slow,  it  is  possible 
to  simplify  the  viscous-unsteady  force  and  avoid  evaluating  the  convolution  integral.  Here  we 
establish  the  condition  under  which  such  a  simplification  can  be  made.  The  main  conclusions 
are:  (1)  Only  when  the  viscous-unsteady  time  scale  (defined  as  the  transition  time  for  the  kernel 
to  change  from  the  slower  decay  to  the  faster  decay)  is  smaller  than  all  the  fluid  time  scales,  the 
viscous-unsteady  formulation  can  be  simplified.  The  resulting  simplified  viscous  unsteady  force 
has  a  form  similar  to  the_added-mass  force.  (2)_In  turbulent  multiphase  flows,  viscous-unsteady 
time  scale  is  larger  than  the  Kolmogorov  time  scale  underjnost  circumstances.  Therefore,  the 
viscous-unsteady  force  generally  needs  to  be  calculated  using  the  convolution  history  integral. 
(3)  For  large  particle-to-fluid  density  ratio,  there  exists  a  narrow  range  of  particle  size  around  the 
Kolmogorov  length  scale,  where  the  simplified  form  of  viscous-unsteady  force  can  be  used. 

Conclusions  on  the  importance  of  unsteady  forces  in  forward  momentum  coupling:  The  ratio 
between  the  magnitudes  of  unsteady  and  quasi-steady  forces  is  used  to  evaluate  the  importance 
of  the  unsteady  forces  in  forward  momentum  coupling.  Key  conclusions  in  regard  to  the 
importance  of  the  unsteady  forces  in  forward  momentum  coupling  include:  (1)  Unsteady  forces 
are  important  in  evaluating  the  motion  of  particles  of  density  comparable  or  smaller  than  that  of 
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the  surrounding  fluid.  (2)  Conventionally,  the  unsteady  forces  are  neglected  in  gas-particle  flows. 
The  present  scaling  analysis  shows  that  for  particles  much  heavier  than  the  fluid  the  importance 
of  the  unsteady  forces  depends  on  the  ratio  between  the  particle  size  and  the  smallest  fluid  length 
scale  and  not  on  the  particle-to-fluid  density  ratio.  (3)  In  DNS  of  turbulent  multiphase  flows,  the 
intent  is  to  accurately  compute  all  the  turbulent  length  scales  and  their  influence  on  particle 
motion.  When  the  particle  size  is  comparable  or  larger  than  the  Kolmogorov  scale  (i.e.,  smallest 
fluid  length  scale),  unsteady  forces  due  to  the  smallest  eddies  are  important  in  accurately 
evaluating  the  particle  motion.  In  LES  and  RANS,  when  the  particle  size  is  larger  than  the 
smallest  resolved  scale,  the  unsteady  forces  due  to  the  resolved  fluid  flow  will  be  important  and 
must  be  included  in  the  computation  of  the  deterministic  component  of  particle  motion.  When 
the  particle  size  is  smaller  than  the  resolved  scale  but  larger  than  the  Kolmogorov  scale,  unsteady 
forces  due  to  the  unresolved  scales  remain  important  and  must  be  accounted  in  the  closure 
models.  (4)  In  shock-particle  interaction,  since  the  particle  size  is  generally  larger  than  the  shock 
thickness  (i.e.,  the  smallest  flow  length  scale),  the  unsteady  forces  dominate  the  quasi-steady 
force  during  shock-particle  interaction.  (5)  For  heavy  particles,  we  observe  the  stress  gradient 
and  added-mass  forces  to  be  of  the  same  order  and  the  ratio  between  viscous-unsteady  and  quasi¬ 
steady  forces  to  scale  as  the  square  root  of  the  ratio  between  the  added-mass  and  quasi-steady 
forces.  Therefore,  in  situations  where  the  inviscid  forces  dominate  the  quasi-steady  force,  they 
also  dominate  the  viscous  unsteady  force. 

Conclusions  on  the  importance  of  quasi-steady  and  unsteady  forces  in  backward  momentum 
coupling:  We  define  momentum  coupling  parameters  as  the  ratio  between  the  back  contribution 
to  the  fluid  momentum  equation  from  the  coupling  force  and  the  fluid  inertial  tenn.  The 
momentum  coupling  parameter  of  each  force  contribution  is  estimated  through  scaling  argument. 
The  key  conclusions  are:  (1)  The  mass  fraction  ratio  is  usually  taken  to  evaluate  the  importance 
of  the  quasi-steady  force  in  backward  momentum  coupling.  It  is  shown  by  the  scaling  analysis 
that  the  momentum  coupling  parameter  of  the  quasi-steady  force  depends  on  the  particle  mass 
and  volume  fraction  ratios  for  heavy  particles  and  bubbles,  respectively.  For  heavy  particles  with 
time  scale  larger  than  the  turbulence  integral  time  scale,  when  Stokes  number  is  greater  than  one, 
the  back  effect  of  quasi-steady  force  on  the  fluid  motion  is  small  and  unsteady  forces  become 
dominant  in  backward  momentum  coupling.  (2)  The  momentum  coupling  parameters  of  stress- 
gradient  and  added-mass  forces  are  similar.  The  contributions  of  the  stress-gradient  and  added- 
mass  forces  to  backward  coupling  are  related  to  the  particle  volume  fraction.  (3)  The  momentum 
coupling  parameter  of  the  viscous-unsteady  force  scales  as  the  square  rOoot  of  the  products  of  the 
momentum  coupling  parameters  of  quasi-steady  and  added-mass  forces.  Therefore,  for  heavy 
particles,  the  back  effect  of  the  viscous-unsteady  force  on  fluid  motion  depends  on  both  particle 
mass  and  volume  fractions. 

4.5  Volume  Fraction  Effect  -  Application  to  Shock-Particle-Curtain  Interaction 
Here  we  built  upon  the  physics-based  interphase  coupling  model  that  has  been  rigorously 
validated  for  an  isolated  particle  subjected  to  shock  waves.  This  work  has  been  performed  in 
collaboration  with  the  Sandia  National  Laboratories.  Experiments  by  Wagner  et  al.  [4], 
conducted  in  the  recently  constructed  Multiphase  Shock  Tube  (MST)  at  the  Sandia  National 
Laboratories,  examined  the  interaction  of  a  planar  shock  wave  with  a  moderately  dense  particle 
curtain.  A  multiphase  shock  tube  was  used  to  drive  a  planar  shock  wave  into  a  particle  curtain, 
i.e.,  a  dense  particle  field  that  was  generated  with  a  gravity-fed  contoured  seeding  method,  which 
is  shown  in  the  schematic  presented  in  Figure  7.  The  initial  volume  fraction  of  particles  in  the 
mixture  was  about  20%.  High-speed  imaging  and  pressure  measurements  showed  that  the 
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interaction  between  the  incident  shock  and  the  particle  curtain  resulted  in  an  upstream- 
propagating  reflected  shock  and  a  downstream-propagating  transmitted  shock.  After  the  passage 
of  the  shock  wave,  the  particle  curtain  was  observed  to  spread  rapidly  as  it  propagated 
downstream.  The  data  suggested  the  expansion  of  the  particle  curtain  to  be  associated  with  a 
difference  in  gas  properties  across  the  streamwise  thickness  of  the  particle  curtain,  although  the 
exact  physical  mechanisms  remained  unexplained  from  the  experiments. 


Fig  7:  Schematic  showing  the  MST  test  section  that  is  used  to  generate  the  dense  particle  curtain.  The 

flow  direction  is  left  to  right. 

The  experimental  conditions  were  simulated  to  investigate  the  underlying  physics  of  shock- 
particle-curtain  interaction  and  to  validate  the  present  compressible  particle-laden  flow  model. 
The  physics  in  the  interaction  between  a  shock  wave  with  a  dense  gas-particle  mixture  is 
markedly  differently  from  that  with  a  dilute  mixture.  Following  the  passage  of  the  shock  wave, 
the  dense  particle  curtain  expands  rapidly  as  it  propagates  downstream.  In  the  simulations,  the 
particles  are  viewed  as  point-particles  and  are  traced  in  a  Lagrangian  framework.  Compared  to 
the  standard  drag  law,  four  major  improvements  are  made  in  the  present  interphase  coupling 
model:  (1)  unsteady  force  contributions  to  particle  force  are  included;  (2)  effect  of 
compressibility  on  hydro-  dynamic  forces  is  included;  (3)  effect  of  particle  volume  fraction  on 
hydrodynamic  forces  is  included;  (4)  effect  of  inter-particle  collision  is  taken  into  account.  The 
complex  behavior  of  the  dense  particle  curtain  is  due  to  the  interplay  between  two-way  coupling, 
finite  particle  inertia,  and  unsteady  forces.  Incorporation  of  these  effects  requires  significant 
modeling  improvements  and  these  improvements  are  observed  to  be  essential  for  the  simulation 
results  to  agree  with  the  experimental  data. 

Figure  9  shows  the  qualitative  comparison  between  Schlieren  images  obtained  from  the 
experiments  and  the  simulations.  More  detailed  quantitative  comparison  is  shown  in  Figure  8, 
where  on  the  x-t  diagram  the  position  of  the  shocks  and  the  particle  curtain  are  compared.  The 
simulations  show  the  rich  details  of  the  particle-laden  flow  arising  from  the  shock-particle- 
curtain  interaction.  The  physical  mechanism  responsible  for  the  expansion  of  the  particle  curtain 
have  been  identified. 
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(a )tf(Lfua)  =  —4.7  ( t  —  —18  fis) 


(b)t/(L/u,)  =  1.3  (t  =  5  fis) 


(c)t/(L/us)  =  5.2  ( t  =  20  fis) 


Upstream  Downstrear 


(d )t/(L/ua)  =  9.1  ( t  =  35  fis) 


{e)t/{L/us)  =  37.2  {t  =  112  fis) 


(f )tf(L/us)  =  70.7  ( t  =  213  fis) 


Fig  8:  Numerically  generated  Schlieren  images  (contours  of  gas  density  gradient  in  region  without 
particles)  are  compared  with  corresponding  experimental  results  for  the  case  of  a  shock  of  Mach  number 
1.92  interacting  with  a  particle  curtain  of  2  mm  thick  at  21%  volume  fraction  of  115  pm  particle.  Results 
are  shown  at  6  different  times.  Close  correspondence  between  the  experiments  and  simulations  can  be 

obsen’ed. 

The  trajectories  of  the  edges  of  the  particle  curtain  computed  with  the  present  model  are  shown 
in  Figure  10  and  are  compared  with  those  computed  with  the  standard  drag  law  of  Clift  and 
Gauvin.37  It  can  be  observed  that  the  upstream  and  downstream  edges  of  the  particle  curtain 
calculated  by  these  two  models  are  substantially  different.  The  present  model  yields  results  that 
match  the  experimental  data  much  better  than  the  standard  drag  law.  If  the  standard  drag  law  is 
used,  both  the  upstream  and  downstream  edges  of  the  particle  curtain  move  slower  than  the 
experimental  results.  For  the  downstream  edge,  at  t  =  150,  the  locations  computed  by  the  present 
model,  the  standard  drag  law,  and  measured  in  the  experiment  are  19.23L,  10.56L,  and  20.47L, 
respectively.  In  other  words,  the  standard  drag  law  underestimates  the  downstream  edge  location 
by  about  50%.  As  a  result  of  the  underestimate  on  the  downstream  edge,  the  expansion  of  the 
particle  curtain  is  also  underestimated  by  the  standard  drag  law.  Therefore,  it  can  be  concluded 
that  if  the  above  improvements  in  the  force  model  are  not  considered,  the  interaction  between  a 
shock  wave  and  a  dense  particle  curtain  and  the  resulting  expansion  of  the  particle  curtain  will 
not  be  captured  accurately.  These  results  are  summarized  in  the  journal  article  J6. 
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Fig  9:  (a)  Comparison  of  the  Reflected  shock  (RS),  Transmitted  shock  (TS),  Upstream  edge  of  particles 
(UEP),  Downstream  edge  of particles  (DEP)  obtained  from  experiments  and  simulations.  Excellent 
agreement  is  seen.  The  importance  of  including  the  different  contributions  to  force  is  demonstrated  in 
frames  (b)  and  (c),  where  the  different  force  contributions  are  plotted  on  short  and  long  time  scales.  Over 
the  initial  short  time  period  (see  frame  b),  during  the  shock  passage  over  the  particle  curtain,  the 
unsteady  forces  ( pressure-gradient ,  added-mass,  viscous  unsteady)  are  much  larger  than  the  quasi-steady 
force.  But  more  interestingly,  even  long  after  the  passage  of  the  shock  the  net  contribution  of  the  unsteady 
forces  remain  comparable  to  the  quasi-steady  contribution.  These  unsteady  forces  are  typically  ignored 
in  standard  modeling  approaches  and  therefore  can  result  in  substantial  error. 


Fig  10:  Trajectories  of  the  upstream  and  downstream  edges  of  particle  curtain  for  Case  1  computed  with 
the  present  model  (“Present")  and  the  standard  drag  law37  (\Std  drag”),  compared  with  the  experimental 

results  (“Exp"). 

4.6  Modeling  and  Simulation  of  Particle-Detonation  Interaction 

Here  we  are  interested  in  the  problem  of  a  strong  shock  interacting  with  a  particle,  where  the 
shock  strength  is  of  sufficient  magnitude  that  the  deformation  and  the  compression  of  the  particle 
cannot  be  ignored.  In  particular,  we  are  interested  in  the  problem  of  detonation-particle 
interaction  in  the  context  of  condensed  phase  detonations,  where  the  density  of  the  products  of 
detonation  is  substantial  and  the  post-shock  pressure  can  be  significantly  larger  than  the  yield 
strength  of  particle.  But  experiments  of  shock  interaction  with  a  deforming  particle  (SIDP)  in  the 
context  of  metal  particles  subjected  to  an  intense  shock  in  a  condensed  matter  ambience  are  rare 
in  the  literature,  due  to  the  complex  nature  of  such  experiments.  Modeling  and  simulation  is  an 
attractive  alternative  that  can  provide  time  resolved  details  of  such  particle  interaction  with  a 
strong  shock  wave.  This  is  the  approach  that  has  been  pursued  here  and  this  work  is  conducted  in 
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collaboration  with  researchers  at  Lawrence  Livermore  National  Laboratory  and  at  the  University 
of  Illinois,  Urbana-Champaign. 

The  primary  goal  of  this  task  is  to  extend  the  earlier  described  theoretical  model  presented  in 
section  4.2  and  propose  a  SIDP  model  that  includes  the  leading  order  effects  of  particle 
deformation  and  compression  heating.  To  reach  this  goal,  we  first  carry  out  DNS  of  a  deforming 
particle  subjected  to  a  strong  shock  in  condensed  matter.  The  material  properties  are 
characterized  by  appropriate  equations  of  state  (e.g.,  Mie-Gruneisen  EOS).  The  DNS  results 
show  the  interaction  details,  such  as  shock  reflection,  refraction  and  evolutions  of  average 
properties  of  particles.  Based  on  the  knowledge  obtained  from  the  DNS  study,  we  improve  the 
model  by  including  the  effects  of  particle  deformation  and  compression  heating.  We  note  that 
shock  interaction  with  a  deforming  particle  (SIDP)  occurs  in  other  interesting  applications  as 
well.  For  example,  shock  interaction  with  a  droplet  or  a  bubble  can  lead  to  significant 
deformation  of  the  droplet  or  the  bubble.  In  fact  a  sufficiently  strong  shock  can  lead  to  violent 
breakup  and  fragmentation  of  the  droplet  or  bubble,  see  e.g.,  [5-7].  The  present  simulation  and 
modeling  effort  was  limited  to  modest  deformation  of  the  particle  and  will  not  include  more 
complex  scenarios  such  as  breakup  and  fragmentation. 


Figure  11  shows  the  DNS  results  of  an  aluminum  particle  subjected  to  an  intense  shock  in 
nitromethane  (shock  speed  is  4800  m/s,  the  post-shock  pressure  and  velocity  are  10.7  Gpa  and 
2270  m/s).  Details  of  particle  deformation,  and  fully  resolved  pressure  and  thermal  fields 
computed  at  the  microscale  are  shown  in  frames  a-d.  The  time  evolution  of  particle  velocity  and 
average  temperature  obtained  from  the  DNS  are  shown  in  Figures  9e  and  9f.  Such  resolution  is, 
however,  impossible  at  the  mesoscale,  where  millions  of  particles  must  be  considered.  So  a 
point-particle  approach  is  employed  at  the  mesoscale,  where  all  the  microscale  details  must  be 
accurately  encapsulated  into  mass,  momentum,  and  energy  coupling  models.  Shown  in  Figure  9e 
and  9f  are  the  particle  velocity  and  temperature  predicted  by  conventional  drag  and  heat  transfer 
point-particle  models  (these  simple  models  are  typically  used  in  meso  and  system-scale 
simulations).  It  is  quite  clear  that  such  standard  models  result  in  very  large  errors,  since  they 
have  not  been  developed  and  tested  for  the  intense  pressure  and  temperature  conditions 
representative  of  detonation  particle  interaction.  As  a  consequence,  it  can  be  expected  that  meso 
and  system-scale  simulations  performed  with  the  standard  models  to  result  in  serious  error.  In 
contrast,  physics-based  coupling  models  developed  by  our  group  can  be  seen  to  perform  quite 
well.  Thus,  Figure  1 1  captures  the  serious  limitations  of  the  standard  micro-to-meso  coupling 
models  and  illustrates  the  importance  of  physics-based  models  rigorously  developed  from 
microscale  simulations  performed  for  extreme  conditions. 


Fig  11:  ALE3D  simulations  of  intense  shock  propagation  over  an  aluminum  particle  embedded  in 
nitromethane  (simulations  performed  by  Dr.  Najjar  at  LLNL).  Frames  (a)  &  (c)  are  pressure  and  (b)  & 
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(cl)  are  temperature  contours  at  0.3  and  0.8  nanoseconds,  respectively.  Frames  (e)  and  (f)  compare  time 
evolution  of  nondimension  al  particle  velocity  and  volume-averaged  temperature  against  standard  and 

improved  physics-based  models. 


Fig  12:  Representative  scenarios  of flow  field  arising  from  shock  interaction  with  a  deformable  particle 

for  different  shock  impedance  and  shock  speed  ratios. 

The  results  of  the  direct  numerical  simulations  have  been  generalized  to  broader  strong-shock- 
particle  interaction  conditions.  Representative  scenarios  of  flow  fields  arising  from  SIDP  are 
shown  in  Figure  12  for  different  combinations  of  the  impedence  ratio  z  /  zm  and  the  shock  speed 

ratio  ust  /  usi  (here  z  ,  zm  are  the  shock  impedence  in  the  particle  and  in  the  surrounding  matrix 

and  ust,  usi  are  the  incidence  and  transmission  shock  propagation  speeds  in  the  matrix  and  in  the 

particle).  The  fonner  ratio  determines  the  reflected  wave  to  be  a  expansion  fan  or  a  shock  wave, 
while  the  latter  ratio  determines  the  shape  of  the  transmitted  shock  (see  Figure  12).  In  summary, 
there  are  four  different  scenarios:  (a)  a  reflected  expansion  fan  and  a  convex  transmitted  shock  ( 
z  ! zm  <  1  and  ust/usi  >1);  (b)  a  reflected  shock  and  a  concave  transmitted  shock  (z  /zm  >  1 

and  ust/usi  <  1);  (c)  a  reflected  expansion  fan  and  a  concave  transmitted  shock  (zp/zm  <  1  and 
ust  / usi  <  1);  and  (d)  a  reflected  shock  and  a  convex  transmitted  shock  (z  / zm  >  1  and  ust  / usi  > 
1).  The  above  presented  are  sample  results.  Further  details  can  be  found  in  the  journal  article  J9. 

4. 7  Instability  of  the  Expanding  Spherical  Gas-Gas  Interface 

Here  we  begin  to  address  the  question  of  how  an  expanding  material  front  becomes  unstable  and 
evolves  into  spatial  patters  seen  in  Figures  1,  2  and  3.  This  instability  and  transition  problem  is  in 
its  infancy.  Rigorous  instability  analysis  of  such  complex  time-dependent  highly-compressible 
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multiphase  flows  has  not  been  performed  in  the  past.  As  the  first  step  we  address  the  instability 
of  a  rapidly  expanding  gas-gas  interface,  before  proceeding  to  the  instability  of  an  expanding 
particulate  front. 

A  high-order  WENO  scheme  is  employed  to  investigate  the  stability  of  a  rapidly  expanding 
material  interface  produced  by  a  spherical  shock  tube.  The  flow  structure  is  characterized  by  a 
forward  moving  primary  shock,  a  backward  moving  secondary  shock,  and  a  spherical  contact 
interface  in-between.  We  consider  herein  the  linear  inviscid  regime  and  focus  on  the 
development  of  the  three  dimensional  perturbations  around  the  contact  interface  by  solving  a 
one-dimensional  system  of  partial  differential  equations.  Numerical  simulations  are  performed  to 
illustrate  the  effects  of  the  contact  interface’s  density  discontinuity  on  the  growth  of  the 
disturbances  for  various  spherical  wave  numbers.  In  a  spherical  shock  tube  the  instability  is 
influenced  by  various  mechanisms  which  include  classical  Rayleigh-Taylor  (RT)  effect,  Bell- 
Plesset  (BP)  or  geometry/curvature  effect,  the  effect  of  impulsively  accelerating  the  interface, 
and  compressibility  effects.  For  an  extended  intermediate  time  period,  it  can  be  shown  that  the 
small  disturbances  grow  exponentially  as  in  the  classical  RT  instability.  During  this  stage,  the 
exponential  growth  rate  increases  with  the  spherical  wave  number,  until  it  saturates  for  very 
large  wave  numbers  due  to  the  finite  thickness  limitation  of  the  numerical  representation  of  the 
contact  interface.  The  results  compare  favorably  with  previous  theoretical  models;  but  indicate 
that  in  addition  to  compressibility,  the  space-time  evolution  of  the  contact  interface’s  thickness 
plays  a  significant  role.  A  parametric  study  is  performed  that  varies  the  pressure  and  density 
ratios  of  the  initial  spherical  container.  The  characteristics  of  the  contact  interface  and  the 
applicability  of  various  instability  theories  is  investigated  for  these  regimes.  Furthermore, 
varying  the  pressure  and  density  ratios  aids  in  understanding  significance  of  compressibility 
effects  on  the  instability  at  different  operating  conditions. 

A  comparison  of  the  computed  growth  rate  with  the  various  theoretical  models  [8-12]  is 
presented  in  Figure  13.  The  predictions  of  planar  (Taylor  [8])  and  spherical  (Plesset  [9],  Espsten 
[12])  RT  instability  theories  are  generally  close  to  each  other  indicating  the  relatively  modest 
effects  of  compressibility  and  sphericity.  The  dominant  feature  of  all  these  theories  is  the 
increase  in  growth  rate  as  square  root  of  the  spherical  wave  number.  The  computed  growth  rate 
follows  this  scaling  up  to  a  wave  numbers  of  about  40,  beyond  which  the  growth  rate  slows 
down  and  eventually  saturates  to  a  constant  growth  rate  that  is  independent  of  the  wave  number. 
In  the  regime  where  growth  rate  is  oc  \pa  the  solution  to  the  compressible  ODE  by  Epstein  [12] 
compare  the  best  with  the  computed  growth  rate.  This  can  be  expected  since  Epstein’s  theory 
accounts  for  both  compressibility  and  radial  effects. 

The  asymptotic  constant  growth  rate  computed  for  large  wave  numbers  compares  well  with 
the  local  convective  growth  rate  given  by  Bandiera  [10].  At  such  large  wave  numbers  the  finite 
thickness  of  the  contact  interface  computed  in  the  base  flow  solution  begins  to  have  a  strong 
influence.  The  interface  does  not  appear  as  a  discontinuity  and  the  density  gradient  in 
conjunction  with  the  local  pressure  gradient  within  the  contact  interface  dictates  the  growth  rate. 
Unlike  classical  RT  instability,  the  convective  instability  is  independent  of  the  spherical  wave 
number.  Figure  9  corresponds  to  a  grid  resolution  of  N  =  16000  for  both  the  base  flow  and 
perturbation  solutions.  The  theory  by  Lindl  [11]  follows  both  the  low  and  large  wave  number 
behaviors.  However,  Lindl’ s  prediction  is  based  on  incompressible  analysis  of  a  planar  interface 
and  is  somewhat  lower  than  the  computed  growth  rate. 
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Fig  13:  The  effect  of  the  wave  number  on  the  growth  rate. 

The  above  inviscid  analysis  has  been  extended  to  the  viscous  regime,  where  the  effect  of  finite 
Reynolds  number  and  thermal  and  mass  diffusivities  on  the  instability  results  have  been 
examined.  Figure  14  demonstrates  the  effect  of  the  Reynolds  number  on  the  growth  rate  for  a 
given  Atwood  number.  For  low  wave  numbers  the  Reynolds  number  has  no  effect  on  the  growth 
rate.  Consequently,  the  growth  rate  predicted  for  wave  numbers  less  than  45  did  not  change  as 
the  Reynolds  number  went  from  6.2  x  105  to  104.  However,  for  high  have  numbers,  as  the 
Reynolds  number  decreases  the  peak  growth  rate  and  its  corresponding  wave  number  also 
decrease.  As  the  Reynolds  number  is  reduced  from  6.2  x  105  to  104,  the  peak  growth  rate  is 
reduced  from  1.15  to  1.07  and  its  corresponding  wave  number  also  reduced  from  400  to  100.  For 
the  limiting  case  of  very  large  Reynolds  number  (Re  — >  oo) ,  the  inviscid  solution  is  recovered 
where  the  growth  rate  saturates  and  does  not  decline  for  very  large  wave  numbers.  The  low- 
Reynolds  number  cases  may  be  taken  as  an  approximation  to  the  turbulent  flow  case,  in  which 
the  presence  of  turbulence  can  be  modeled  by  an  eddy  viscosity  that  increases  the  effective 
Reynolds  number.  Further  details  can  be  found  in  the  articles  J4  and  J8. 
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Fig  14:  The  effect  of  the  Reynolds  number  on  the  growth  rate  for  an  Atwood  number  of  0.3575. 
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4.8  Instability  of  a  Rapidly  Expanding  Spherical  Particulate  Interface 

The  problem  under  study  is  that  of  a  mixture  of  gas  and  metal  particles  (e.g.  aluminum)  initially 
contained  in  a  spherical  diaphragm  of  a  small  radius  under  high  pressure.  The  gas-particle 
mixture  is  then  suddenly  released  as  in  a  spherical  shock-tube.  The  problem  represents  several 
real-life  applications,  such  as,  particle-laden  detonative  combustion,  dispersal  of  metalized 
explosives,  volcanic  eruption,  and  ignition  of  solid-propellant  rockets,  to  mention  only  a  few.  In 
the  single-phase  case,  when  the  high-  pressure  gas  is  suddenly  released  from  a  spherical 
container,  a  forward-moving  shock  wave  is  formed  along  with  a  backward-moving  rarefaction 
wave,  whose  trailing  edge  becomes  a  secondary  shock  due  to  the  convergent  geometry.  Between 
the  two  shocks  a  contact  interface  (Cl)  is  formed,  which  acts  as  a  material  interface  that  separates 
the  hot  driver  gas  from  the  cold  driven  pure  gas.  Due  to  the  inherent  density  jump  across  this 
contact  interface,  Rayleigh-Taylor-Instability  (RTI)  emerges  where  inter-penetrating  dimples 
form,  around  which  there  is  considerable  vorticity.  In  the  case  of  a  gas-particle  mixture  in  a 
spherical  shock-tube,  a  second  interface  is  fonned  that  separates  the  gas-particle  mixture  from 
the  pure-gas  phase.  The  presence  of  the  solid  particles  influences  the  instability  of  the  gas-gas 
Cl,  which  is  the  subject  of  the  present  study. 

In  the  multiphase  context,  the  effect  of  the  presence  of  metal  particles  on  the  gas-gas  Cl  needs  to 
be  studied.  The  gas-particle  interactions  are  quite  complicated  because  there  are  several 
additional  mechanisms  involved,  such  as  the  quasi-steady  viscous  forces,  the  heat  transfer 
between  the  particles  and  the  gas,  and  the  added-mass  forces.  The  instability  of  the  gas-gas  Cl  is 
affected  by  the  presence  of  the  particles  at  two  levels.  In  the  first,  the  presence  of  the  particles 
modifies  the  gas  base  flow,  which  drives  the  Cl  instability.  At  the  second  level,  the  gas-particle 
interaction  produces  additional  perturbation  source-tenns  in  the  stability  equations  of  the 
perturbed  gas  flow.  Although  an  Eulerian-  Lagrangian  approach,  if  utilized  properly,  can  capture 
the  effect  of  particles  on  the  contact  interface,  it  is  harder  to  capture  the  instability  of  the  particle 
interface  itself  if  needed.  Hence,  a  two-fluid  Eulerian  approach  is  implemented  herein.  The 
methodology  employed  in  Mankbadi  &  Balachandar  (2012)  for  the  single-phase  problem  is 
extended  here  to  investigate  the  complex  multiphase  problem. 

The  base  flow’s  discontinuous  surfaces  are  tracked  in  the  presence  of  solid  aluminum  particles. 
The  density  of  aluminum  particles  relative  to  that  of  the  initial  compressed  gas  yields  a  ratio  of 
43.66.  In  Figure  15,  the  locations  of  the  interfaces  for  three  different  combinations  of  particle 
diameter  and  volume  fraction  are  depicted.  Figure  2  (A)  (dp  =  10  micron,  initial  particle  volume 
fraction  =  1.5%)  illustrates  the  locus  of  the  primary  shock  (PS),  the  secondary  shock  (SS),  the 
contact  interface  (Cl)  separating  the  driver  gas  from  the  driven  gas,  and  the  particles  interface 
(PI).  In  the  near  field,  the  secondary  shock  initially  moves  out  but  soon  turns  and  travels  inward. 
For  very  early  time  the  secondary  shock  is  in  close  proximity  to  both  the  contact  interface  and 
the  particle  interface.  The  particle  interface  impacts  the  behavior  near  the  contact  interface, 
especially  when  it  collides  with  the  contact  interface  and  overtakes  it  at  about  t=2.5.  The 
particles  then  pierce  through  the  contact  interface  and  are  no  longer  in  the  initially  compressed 
gas  but  rather  are  now  in  the  shock-heated  fluid.  Now,  consider  two  other  extreme  cases  for  the 
base  flow.  In  Figure  15  (b)  ( dp  =  500  micron,  initial  particle  volume  fraction  =  5%),  the  particle 
response  time  has  become  too  large,  and  as  such  the  particles  are  not  responding  quickly  to  the 
gas  flow.  In  this  case,  the  PI  stays  away  from  the  Cl  as  shown  in  the  diagram.  In  the  other 
extreme  (dp  =  1  micron,  initial  particle  volume  fraction  =  5%),  the  particle  response  time  has 
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become  very  small  such  that  particles  respond  quickly  to  the  gas,  and  consequently  the  PI  and  Cl 
are  close  to  each  other,  as  illustrated  in  Figure  15  (c). 


*  (a) 


6  (b)  5  (c) 

Fig  15:  The  T  diagram  depicting  the  radial  location  of  the  primary  shock,  particle  interface,  contact 
interface,  and  secondary  shock.  A)  For  dp  —  10  pm,  pp/Pg  —  43.66,  and  0P()  =  0.015.  B)  For 
dp  —  500pm,  pp/pg  —  43.66,  and  <pPo  =  0.05.  C)  For  dp  —  1pm,  pp/pg  —  43.66,  and  <pPo  = 

0.05. 

4.9  Non-linear  Instability  of  the  Expanding  Cylindrical  Gas-Gas  Interface 
When  the  amplitude  of  the  perturbation  grows  to  10-20%  of  the  wavelength,  non-linear  effect 
becomes  important.  The  symmetry  between  the  light  and  heavy  fluids  breaks  down  and  a  bubble- 
spike  structure  appears.  Conventionally,  bubble  refers  to  the  light  fluid  that  rises  into  the  heavy 
fluid,  while  spike  denotes  the  heavy  fluid  that  penetrates  into  the  light  fluid.  The  non-linear 
growth  is  characterized  by  the  heights  of  the  bubble  and  spike.  In  the  non-linear  regime,  it  has 
been  shown  previously  using  simulations,  theoretical  analysis  and  experiments  that  the  heights  of 
the  bubbles/spikes  increase  at  a  constant  velocity  when  a  single  mode  perturbation  is  introduced 
at  the  interface.  In  other  words,  the  perturbation  amplitude  grows  linearly  in  time.  Recently, 
investigations  show  that  the  growth  velocities  of  the  heights  of  the  bubbles/spikes  increase  again 
after  being  at  a  constant  value  for  a  finite  period  of  time. 

Modeling  of  the  non-linear  growth  of  RTI  has  been  intensively  studied  in  the  past.  The  goals  of 
the  present  study  are  to  improve  the  understanding  of  RTI  in  explosion  flows  through  rigorous 
numerical  simulations  of  the  present  problem  and  to  examine  the  capabilities  of  the  existing  RTI 
models  on  predicting  RTI  in  explosion  flows.  A  simple  extension  to  Miles’  model  is  proposed,  in 
which  the  effects  of  the  transition  from  linear  to  non-linear  growth  is  captured.  In  addition,  the 
background  flow  density  variations  are  included.  As  can  be  seen  from  Figure  16(a)  the  RTI  can 
still  be  studied  as  bubbles  and  spikes  even  at  late  times.  The  mushroom- type  structures  are  still 
intact.  So  we  can  conclude  that,  for  a  given  wavenumber,  the  non-linear  growth  in  the  azimuthal 
and  longitudinal  directions  are  quite  different.  While  in  the  azimuthal  case,  the  structures  are 
well  behaved  for  large  wavenumbers  even  for  late  times,  severe  mixing  ensues  for  the 
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corresponding  longitudinal  case  (see,  Figure  16(b)).  This  shows  that  for  a  given  wavenumber,  the 
non-linear  growth  does  depend  on  the  direction  in  which  the  mode  is  imposed. 

y 

V 

Fig  16:  (a)  Bubble  and  spike  mushroom-type  structures  for  azimuthal  wavenumber  (m  =  48) 
at  time  =  6:9  and  (b)  No  clear  bubble-spike  pattern  for  longitudinal  wavenumber  (n  =  48)  at  t  =  6:9 

5.  Recommendations  and  Future  Work 

We  are  continuing  the  research  along  all  the  fronts  outlined  in  section  4.  In  particular,  we  plan  to 
accomplish  further  fundamental  advances  in  interphase  coupling  modeling,  experimental 
validation,  and  prediction  of  interfacial  instability  of  explosive  multiphase  material  fronts.  These 
ongoing  activities  are  discussed  below: 

Modeling:  Despite  the  good  performance  of  the  present  interphase  coupling  model  described 
above,  there  are  still  several  complex  scenarios  of  practical  importance,  for  which  the  coupling 
model  requires  further  investigation  and  improvement. 

1)  In  the  limit  of  moderate  particle  clustering,  the  effect  of  local  particle  volume  fraction  can  be 
expected  to  play  a  key  role  in  the  quasi-steady,  and  the  unsteady  force  contributions.  Currently 
we  rely  upon  incompressible  flow  results  to  provide  up  guidance  for  modeling  the  volume 
fraction  effect  in  the  compressible  flow  regime.  We  will  pursue  direct  numerical  simulations  at 
the  microscale  of  steady  and  unsteady  flows  around  periodic  array  of  particles  to  develop 
fundamental  understanding  of  volume  fraction  effect  on  quasi-steady  and  unsteady  forces  in  the 
compressible  flow  regime. 

2)  When  particle  volume  fraction  is  large,  then  particle-particle  interaction  in  the  fonn  of  inter¬ 
particle  collisions  will  begin  to  play  an  important  role.  In  [J6]  we  used  a  simple  particle-pressure- 
based  model  to  account  for  the  inter-particle  collision  force.  This  simple  representation  was 
sufficient  for  the  modest  volume  fraction  considered  in  that  investigation.  In  order  to  extend  the 
formulation  to  higher  volume  fraction,  more  rigorous  understanding  of  inter-particle  collision  is 
required,  especially  in  the  context  of  a  shock  wave  propagating  over  a  dense  curtain  of  particles. 
On  one  hand  we  are  coordinating  with  researchers  at  Sandia  National  Laboratories  to  extend 
their  experiments  [4]  to  denser  particle  curtains,  where  inter-particle  collision  will  be  the 
dominant  force  on  the  particles.  A  rigorous  universal  compressible  multiphase  formulation  is 
currently  under  development.  We  have  started  to  first  explore  its  mathematical  structure  and 
limiting  analytical  and  numerical  behaviors.  As  we  complete  this  theoretical  framework,  we  are 
continuing  to  test  the  formulation  by  comparing  its  predictions  against  high  quality  experimental 
results,  including  those  from  the  multiphase  shock  tube  facility  at  Sandia  National  Laboratories. 
We  are  begining  to  perform  Discrete  Element  Method  (DEM)  simulations  along  with  our  point- 
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particle  coupling  model  of  the  hydrodynamic  forces.  The  DEM  results  can  be  used  to  develop 
improved  inter-particle  collision  force  models  and  to  validate  them  against  experimental 
measurements.  However,  in  the  DEM  approach  an  empirical  model  is  required  in  accounting  for 
the  momentum  and  energy  exchange  during  particle -particle  contact. 

3)  Note  that  once  the  detonation  wave  completely  sweeps  across  the  entire  explosive,  the 
detonation  wave  will  be  transmitted  as  a  shock  wave  (blast  wave)  through  the  surrounding  air.  At 
the  initial  stages  the  intensity  of  this  shock  wave  is  very  high  and  the  corresponding  shock  Mach 
number  can  be  order  ten.  The  compressible  coupling  models  that  have  been  developed  and  tested 
have  generally  been  limited  to  subsonic  and  transonic  conditions.  Thus,  there  is  a  strong  need  to 
evaluate  the  appropriateness  of  the  interphase  coupling  models  at  even  extreme  flow  conditions, 
and  develop  new  models  (or  modifications)  as  necessary. 

4)  We  are  continuing  our  investigation  of  detonation-particle  interaction,  where  the  particle 
undergoes  significant  deformation  (see  discussion  in  section  4.6).  We  are  performing  ALE3D 
simulations  of  detonation-particle  interaction  under  varying  conditions  of  incident  shock/ 
detonation  wave  intensity,  particle-to-fluid  density  and  impedance  ratios.  These  simulation 
results  are  being  used  to  develop  a  more  robust  model  that  fully  account  for  particle  deformation 
into  consideration.  Efforts  are  underway  to  extend  the  point-particle  models  to  predict  not  just 
the  mass-averaged  velocity  and  heating  but  to  predict  the  peak  values  of  stress  and  temperature 
as  well. 

Experimental  Validation:  Rigorous  validation  is  essential  at  every  stage  of  model  development 
and  integrated  simulation.  We  have  taken  care  to  find  appropriate  validation  data  for  each 
component  of  the  coupling  model,  as  well  as  their  overall  integrated  contribution.  As  we 
incorporate  more  features  into  the  model  additional  experimental  results  are  being  identified  for 
rigorous  validation. 

1)  The  present  very-strong-shock  or  detonation-particle  interaction  model  has  been  motivated  by 
detailed  numerical  simulations  performed  at  the  microscale.  However,  the  model  must  be 
validated  against  high  quality  experimental  results.  Measurements  of  deformation  and  motion  of 
an  isolated  particle  resulting  from  its  interaction  with  a  detonation  wave  is  a  challenging  task. 
We  are  currently  collaborating  with  researchers  at  the  Air  Force  Research  laboratory  -  Munitions 
Directorate  to  perform  this  kind  of  experiments  along  with  detailed  time-resolved  x-ray  and  other 
high-speed  measurements. 

2)  Prof.  D.  Frost  at  McGill  University  and  his  collaborators  have  performed  outstanding 
collection  of  experiments  on  explosive  dispersal  of  particle  under  varying  conditions.  Prof.  Frost 
is  current  on  his  sabbatical  visit  at  our  group  at  the  University  of  Florida.  We  are  closely 
collaborating  with  him  and  applying  our  modeling  and  simulation  approach  to  his  experimental 
configurations. 

3)  Interphase  energy  coupling  models  of  shock-particle  and  detonation-particle  interaction  have 
hardly  been  validated  due  to  paucity  of  high-quality  experimental  results.  Thus,  validation  of 
thennal  coupling  models  is  of  great  interest.  In  particular,  we  desire  validation-quality 
experimental  measurements  of  particle  and  gas  temperature  immediately  following  a  strong- 
shock  or  detonation  particle  interaction. 

Instability:  Our  investigation  of  instability  of  a  rapidly  expanding  material  interface  presented  in 
section  4.5  was  limited  to  only  gas-gas  interface.  Thus  we  only  considered  the  case  of  a  high- 
density  highly  compressed  sphere  of  gas  (for  example  products  of  a  spherical  detonation) 
expanding  into  the  ambient.  Due  to  the  near-constant  deceleration  of  the  gas-gas  interface  it 
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initially  undergoing  Rayleigh-Taylor  (RT)  instability,  which  was  investigated  both  in  the  inviscid 
and  viscous  limits  (see  section  4.5).  Our  analysis  has  established  a  firm  mathematical 
methodology  for  investigating  instability  of  such  rapidly  expanding  material  fronts.  Over  the 
next  two  years  we  plan  to  exploit  this  methodology  and  investigate  the  following  extensions. 

1)  We  are  extending  the  investigation  of  the  instability  of  a  rapidly  expanding  particulate  front, 
presented  in  section  4.8.  We  will  also  consider  cases  similar  to  that  in  Figure  2e,  where  the 
particles  were  placed  around  an  explosive.  The  introduction  of  the  particulate  phase  in  the 
instability  analysis  introduces  additional  complexities,  and  the  role  of  interphase  momentum  and 
energy  coupling  will  be  of  importance. 

2)  One  of  the  interesting  features  of  the  spherical  instability  is  that  the  growth  rate  only  depends 
on  the  spherical  model  number  n  and  does  not  distinguish  between  the  combination  of 
circumferential  and  azimuthal  components.  We  plan  to  continue  performing  the  non-linear 
analysis  of  the  instability  through  direct  numerical  simulations,  under  cylindrical  configuration 
(similar  those  presented  in  section  4.9). 
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